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In  this  study,  a  general  model  of  proton  exchange  membrane  fuel  cell  (PEMFC)  was  constructed,  imple¬ 
mented  and  employed  to  simulate  the  fluid  flow,  heat  transfer,  species  transport,  electrochemical  reaction, 
and  current  density  distribution,  especially  focusing  on  liquid  water  effects  on  PEMFC  performance.  The 
model  is  a  three-dimensional  and  unsteady  one  with  detailed  thermo-electrochemistry,  multi-species, 
and  two-phase  interaction  with  explicit  gas-liquid  interface  tracking  by  using  the  volume-of-fluid  (VOF) 
method.  The  general  model  was  implemented  into  the  commercial  computational  fluid  dynamics  (CFD) 
software  package  FLUENT®  v6.2,  with  its  user-defined  functions  (UDFs).  A  complete  PEMFC  was  consid¬ 
ered,  including  membrane,  gas  diffusion  layers  (GDLs),  catalyst  layers,  gas  flow  channels,  and  current 
collectors.  The  effects  of  liquid  water  on  PEMFC  with  serpentine  channels  were  investigated.  The  results 
showed  that  this  general  model  of  PEMFC  can  be  a  very  useful  tool  for  the  optimization  of  practical 
engineering  designs  of  PEMFC. 

©  2008  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

In  recent  years,  high-performance  computing  and  advanced 
numerical  algorithms  have  allowed  researchers  to  model  proton 
exchange  membrane  fuel  cell  (PEMFC)  systems  as  well  as  individual 
components  with  greater  fidelity  than  ever  before.  In  gen¬ 
eral,  PEMFC  operations  involve  simultaneously  multi-component, 
multi-phase,  multi-dimensional  fluid  flow  with  heat  and  mass 
transfer  and  electrochemical  reactions.  Hence,  a  complete  math¬ 
ematical  model  is  necessary  to  characterize  more  fully  the  physical 
behavior,  to  aid  our  understanding  of  complex  phenomena  occur¬ 
ring  in  a  fuel  cell  system  and  to  provide  powerful  tools  for  fuel  cell 
design  and  optimization. 

Fuel  cell  performance  is  evaluated  by  the  relation  between  the 
cell  voltage  and  current  density.  The  voltage  loss  is  mainly  caused 
by  the  activation  loss  due  to  electrochemical  reactions,  the  Ohmic 
loss  due  to  resistance  and  fuel  crossover  and  the  mass  transport 
loss  due  to  limitation  of  gas  transport  inside  the  cell.  During  fuel 
cell  operation,  by-product  water  is  generated  in  the  cathode  cata¬ 
lyst  layer  in  the  form  of  liquid,  leading  to  a  gas-liquid  flow  in  the 
porous  media  and  the  channels.  On  the  other  hand,  due  to  the  low 
operating  temperature  of  PEMFCs  (30-100  °C),  excessive  humid¬ 
ification  could  result  in  water  vapor  condensation.  These  could 
subsequently  block  the  gas  flow  channels  resulting  in  a  lower  air¬ 
flow  rate  on  the  cathode  side,  increasing  the  voltage  loss  due  to 
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mass  transport,  thus  decreasing  fuel  cell  performance.  However, 
due  to  the  special  chemical  structure  of  the  proton  exchange  mem¬ 
brane  (PEM),  the  membrane  must  be  well  hydrated  to  ensure  that 
a  sufficient  amount  of  hydrogen  ions  could  cross.  Water  content 
is  also  an  important  factor  that  affects  the  Ohmic  resistance  in 
the  membrane.  Therefore,  keeping  an  appropriate  amount  of  water 
content  in  the  fuel  cell  to  avoid  both  membrane  dehydration  and 
water  vapor  condensation  has  been  a  critical  issue  in  improving 
fuel  cell  performance.  In  reality,  however,  it  is  almost  impossible 
to  manage  water  on  both  the  anode  and  cathode  sides  without 
dehydration  and  condensation;  this  is  simply  because  water  vapor 
condensation  in  the  gas  flow  channels  of  practical  fuel  cell  appli¬ 
cations  is  unavoidable.  Therefore,  water  management,  especially 
liquid  water  management,  to  which  many  engineers  and  scientists 
have  recently  paid  particular  attention,  has  been  a  critical  challenge 
for  a  high-performance  fuel  cell  design  and  optimization. 

In  the  last  decade,  water  management-related  studies  were  per¬ 
formed  numerically  and  experimentally  for  different  purposes  and 
in  several  ways.  Recently,  three-dimensional  models  based  on  the 
computational  fluid  dynamics  (CFD)  approach  have  been  presented 
by  different  commercial  CFD  software  packages  such  as  FLUENT®, 
CFX®,  Star-CD®,  and  CFDRC®.  A  CFD  modeling  of  PEMFCs  which 
simultaneously  considered  the  electrochemical  kinetics,  current 
distributions,  hydrodynamics,  and  multi-component  transport  was 
conducted  by  Um  et  al.  [1].  A  three-dimensional  numerical  simu¬ 
lation  of  a  straight  gas  flow  channel  in  a  PEMFC  was  performed  by 
Dutta  et  al.  [2]  using  a  commercial  CFD  software  FLUENT®.  Hon- 
tanon  et  al.  [3]  also  employed  FLUENT®  to  implement  their  3D, 
stationary  gas  flow  model.  A  study  exploring  the  steady-state  gas 
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Nomenclature 

a 

water  activity 

As 

heat  transfer  surface  area,  m2 

ASUrf 

reactive  surface  area,  m2 

Cf 

concentration  of  sulfonic  acid  ions 

(HS03-), 

kmol  m3 

CP 

specific  heat  capacity,  J  kg-1  K-1 

Q 

species  concentration  i,  kmol  m-3 

Di 

diffusion  coefficient  of  species  i  in  gas 

mixture, 

m2  s-1 

F 

Faraday  constant,  9.6487  x  107  Ckmol-1 

h 

convective  heat  transfer,  W  m-2  K-1 

hH20 

the  enthalpy  of  formation  of  water  vapor,  N  m  kg-1 

I 

current  density,  Am-2 

Gve 

average  current  density,  Am-2 

J 

mass  flux,  kg  m-2  s-1 

ke  ff 

effective  thermal  conductivity,  W  m-1  K-1 

Mi 

molecular  weight  of  species  i  in  gas 

mixture, 

kg  kmol-1 

nd 

electro-osmotic  drag  coefficient 

nf 

charge  number  of  the  sulfonic  acid  ion 

P 

pressure,  Pa 

Pi 

partial  pressure  of  species  i,  Pa 

R 

universal  gas  constant,  8314Jkmol-1  K-1 

Q 

heat  rate,  W 

Rc at,  Ran 

volumetric  current  density,  Am-3 

S 

source  term 

t 

time,  s 

T 

temperature,  K 

U,  V,  w 

velocities  inX,  Y,  and  Z  directions,  respectively,  m  s-1 

Voc 

open-circuit  potential,  V 

Veen 

cell  potential,  V 

Vref 

reference  potential,  V 

4^ 

volume,  m3 

Xf 

mole  fraction  of  species  i 

Yi 

mass  fraction  of  species  i 

Greek  symbols 

a 

transfer  coefficient 

P 

the  factor  accounts  for  energy  release 

y 

concentration  dependence 

Yp'Yt 

exponent  factors 

£ 

porosity 

0W 

contact  angle,  ° 

r\ 

overpotential,  V 

K 

surface  curvature 

hydraulic  permeability,  m2 

K(j) 

electrokinetic  permeability,  m2 

X 

water  content 

ll 

dynamic  viscosity,  kgm-2  s-1 

p 

density  of  gas  mixture,  kgm-3 

Pi 

density  of  species  i,  kg  m-3 

o 

phase  conductivity,  m-1 

X 

gaseous  permeability,  m2 

V 

reaction  rate,  kmol  m2  s-1 

0 

phase  potential,  V 

<p 

relative  water  content 

X 

surface  tension  coefficient,  N  m-1 

CO 

excess  coefficient 

Subscripts  and  superscripts 

an 

anode 

cat 

cathode 

e  electrochemical  reaction 

eff  effective 

g  gas  phase 

H2  hydrogen 

H20  water 

i  species  i 

in  inlet 

1  liquid  phase 

m  membrane  phase 

out  outlet 

N2  nitrogen 

02  oxygen 

ref  reference 

s  solid  phase 

sat  saturated 

surf  surface 

w  water  vapor 


transport  phenomena  in  micro-scale  parallel  flow  channels  was 
conducted  by  Cha  et  al.  [4]  in  which  oxygen  concentration  along 
a  single  gas  flow  channel  and  other  flow  patterns  that  may  affect 
fuel  cell  performance  were  discussed.  Similarly,  gas  concentration 
of  a  steady-state  flow  along  fuel  cell  flow  channels  was  obtained 
numerically  by  Kulikovs ky  [5].  However,  in  all  the  studies  men¬ 
tioned  above,  the  effects  of  liquid  water  were  neglected.  Yi  et  al.  [6] 
pointed  out  that  water  vapor  condensation  was  inevitable  on  both 
the  anode  and  cathode  sides  of  a  PEMFC,  and  they  discussed  a  liquid 
water  removal  technique  that  used  a  water  transport  plate  to  lead 
excess  liquid  water  to  the  coolant  flow  channels  by  a  pressure  dif¬ 
ference.  Wang  et  al.  [7]  conducted  a  two-phase  model  on  PEMFC 
cathode  to  address  the  liquid  water  saturation.  You  and  Liu  [8] 
also  considered  liquid  water  saturation  in  a  straight  channel  on  the 
cathode  side.  Both  the  Refs.  [7,8]  showed  the  importance  for  consid¬ 
ering  liquid  water  in  numerical  modeling  of  PEMFCs.  A  3D  model 
by  Natarajan  and  Nguyen  [9]  introduced  the  effects  of  flow  rate, 
inlet  humidity  and  temperature  on  the  liquid  water  flooding  in  the 
PEMFC  cathode.  The  numerical  results  of  this  model  were  validated 
with  experimental  data.  Berning  and  Djilali  [10]  presented  a  multi¬ 
phase,  multi-component  3D  model  with  heat  and  mass  transfer, 
where  liquid  water  transport  in  GDLs  was  numerically  modeled  by 
using  viscous  and  capillary  effects.  This  method  was  also  used  in 
a  3D  model  by  Mazumder  and  Cole  [11].  Other  approach  to  deal 
with  two-phase  flow  in  GDLs  based  on  thermodynamic  equilib¬ 
rium  conditions  was  given  by  Vynnycky  and  Vynnycky  [12].  In  this 
approach,  the  location  of  interface  between  one  phase  and  the  other 
phase  was  pointed  out  from  its  numerical  results.  In  recent  years, 
more  two-phase  models  have  been  published  [13-15],  these  sim¬ 
ulations  predicted  water  flooding  inside  PEMFCs,  and  the  liquid 
water  effects  on  PEMFC  performances.  Large-scale  simulations  for 
complex  flow  field  were  also  performed  with  experimental  valida¬ 
tions  [15-18]. 

All  the  above  research  papers  have  significantly  contributed  to 
PEMFC  research  and  development.  However,  by  far,  to  the  authors’ 
knowledge,  most  of  the  two-phase  numerical  models  have  not  con¬ 
sidered  the  interface  tracking  between  liquid  water  and  gas.  The 
detailed  behaviors  of  liquid  water  transport  inside  the  PEMFCs  were 
rarely  discussed  except  for  the  present  authors’  previous  study  [19], 
which  only  dealt  with  part  of  serpentine  channels— the  single  U- 
shaped  channel.  Recently,  Zhou  et  al.  [20-22]  also  conducted  two 
more  studies  that  dealt  with  liquid  water  in  serpentine  and  straight 
parallel  fuel  cell  stacks.  The  results  showed  that  different  designs  of 
gas  diffusion  layers  (GDLs)  and  flow  channels  will  affect  the  liquid 
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water  flow  patterns  significantly,  thus  influencing  the  performance 
of  PEMFCs.  Djilali  and  co-workers  [23]  presented  a  2D,  numerical 
investigation  of  the  dynamic  behavior  of  liquid  water  entering  a 
PEMFC  channel  through  a  small  hole  that  was  assumed  as  a  GDL 
pore.  In  these  studies,  the  water  behaviors  can  give  very  useful 
insights  about  water  management  only  in  the  cathode  channels 
although  the  electrochemical  reactions  were  not  considered. 

Numerical  models  that  considered  liquid  water  in  the  porous 
mediums  such  as  GDL  were  developed  in  several  ways.  Nam  and 
Kaviany  [24]  developed  a  two-dimensional,  two-phase  numeri¬ 
cal  model  by  considering  random  carbon  fiber  mats  as  the  GDL. 
Single-  and  two-layer  diffusion  mediums  were  both  considered  to 
investigate  the  effective  diffusivity  and  water  saturation.  Two  main 
factors  influencing  on  two-phase  transport  in  GDL  of  a  PEMFC  such 
as  permeability  and  capillary  pressure  were  also  determined  via 
a  capillary  network  model  developed  by  Djilali  and  co-workers 
[25].  In  addition,  the  GDL  considered  as  a  network  of  pore  bod¬ 
ies  and  pore  throats  used  to  address  the  liquid  water  behavior  in 
PEMFC  were  presented  by  Fowler  and  co-workers  [26]  or  Sinha 
and  Wang  [27].  A  study  on  the  interaction  between  the  GDL  and 
the  flow  field  was  performed  by  Dohle  et  al.  [28]  numerically  and 
experimentally.  Other  models  that  considered  the  porous  medi¬ 
ums  also  mainly  focused  on  the  porosity  of  the  carbon  fiber  paper 
that  could  influence  the  performances  of  PEMFCs  [29,30].  How¬ 
ever,  the  detailed  flow  patterns  that  liquid  water  exhibits  across 
the  porous  medium  of  GDL  was  rarely  discussed.  Based  on  their 
previous  efforts  on  water  behaviors  inside  cathode  channels,  the 
present  author’s  group  further  investigated  water  behaviors  inside 
innovative  GDLs  [31,32]. 

Based  on  the  above  literature  reviews,  it  is  evident  that  there 
is  an  urgent  need  to  develop  a  general  mathematical  model 
with  all  detailed  physics  included,  e.g.,  multi-phase  with  VOF 
interface-tracking  between  gas-phase  and  liquid-water-phase, 
multi-components,  heat  and  mass  transfer,  electrochemical  reac¬ 
tions,  and  water-phase-change  effects.  In  this  paper,  a  general 
model  is  presented  with  all  parts  of  PEMFC  including  membrane, 
catalyst  layers,  GDLs,  gas  flow  channels  and  current  collectors. 

2.  Mathematical  model 

Fig.  1  schematically  shows  a  disassembling  diagram  of  a  PEMFC 
comprising  of  two  distinct  current  collectors  with  flow  channels 
formed  on  them  at  both  the  cathode  and  anode;  an  MEA  consists 
of  a  proton  exchange  membrane  sandwiched  between  two  catalyst 
layers.  Between  the  current  collectors  and  the  MEA,  there  are  two 
GDLs  on  both  sides.  For  the  sake  of  simplicity  and  popularity,  the 
flow  channels  are  in  serpentine  form.  In  this  paper,  the  mathemati¬ 
cal  model  considers  all  the  parts  shown  in  Fig.  la,  the  computation 
domain  is  shown  in  Fig.  lb,  with  geometrical  parameters  listed  in 
Table  1 .  The  computation  domain  consists  of  current  collectors,  flow 
channels,  GDLs,  catalyst  layers  and  membrane.  The  channels  are 
serpentine,  total  path  length  for  a  channel  of  115  mm  with  1  mm 
(width)  x  1  mm  (height)  in  the  cross-section  flow  area.  The  species 
considered  are  hydrogen,  oxygen,  water  vapor  and  nitrogen. 

2.1.  Model  assumptions 

The  assumptions  used  in  developing  the  model  are  as  follows: 

1.  Ideal  gas  law  was  employed  for  gaseous  species. 

2.  The  fluid  flow  in  the  fuel  cell  was  laminar  due  to  the  low  flow 

velocities  and  the  small  size  of  gas  flow  channels. 

3.  The  porous  media  including  membrane,  catalyst  layers  and  GDLs 

were  considered  to  be  isotropic. 


Fig.  1.  (a)  Schematic  diagram  of  PEMFC  and  (b)  computation  domain. 


Table  1 

Physical  properties  and  parameters 


Physical  properties  or  parameters 

Value 

Total  channel  length 

0.115  m 

Channel  width 

0.001  m 

Channel  height 

0.001  m 

Membrane  thickness 

50  x 10~6  m 

GDL  thickness 

300  x 10~6  m 

Catalyst  layer  thickness 

10  x 10-6  m 

Anode  inlet  pressure,  Pan  in 

2  atm 

Cathode  inlet  pressure,  Pcatijn 

2  atm 

Anode  and  cathode  inlet  temperature 

300  K 

Anode  inlet  excess  coefficient,  coa 

3 

Cathode  inlet  excess  coefficient,  coc 

3 

Relative  humidity  of  anode  inlet 

20% 

Relative  humidity  of  cathode  inlet 

100% 

Open-circuit  voltage,  Voc 

1.15  V 

Faraday  constant,  F 

96,487,000  Ckmol”1 

Gas  constant,  R 

8314Jkmol_1  K-1 

Anode  volumetric  reference  exchange  current 

7  x  1010  Akmol-1 

density/reference  hydrogen  concentration, 

Cathode  volumetric  reference  exchange  current 

7  x  105  Akmol-1 

density/reference  oxygen  concentration, 

«£/«%)** 

Anode  transfer  coefficient,  o;an 

0.5 

Cathode  transfer  coefficient,  aCat 

0.5 

Anode  concentration  dependence,  yan 

0.5 

Cathode  concentration  dependence,  yc at 

1.0 

Factor  accounts  for  energy  release,  ft 

0.5 

Membrane  porosity,  £mem 

0.5 

Diffusion  layer  porosity,  £Gdl 

0.5 

Catalyst  layer  porosity,  £cataiySt 

0.5 

Permeability  of  porous  media,  r 

1.76  x  10-11  m2 

Contact  angle,  0W 

90° 

Surface  tension,  x 

0.065  Nm-1 2 3 
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4.  The  fuel  cell  cooling  was  controlled  by  forced  convection  heat 
transfer. 

2.2.  Governing  equations 

22  A.  Mass  conservation  equation 

The  continuity  equation  is  expressed  as  follows: 

^  +  V-(spv)  =  Sm  (1) 

For  multi-phase  flow,  the  mixture  density  is  defined  as  follows 
[33,34]: 


P  =  S\p\  +  (1  —  5\)pg 
s\  +  sg  —  1 


(2) 

(3) 


2.2.2.  Momentum  conservation  equation 

The  momentum  equation  is  dependent  on  the  volume  fractions 
of  all  phases: 

9 

^(gpv)  +  V(£pvv)  =  -SVp  +  V[£/xVv]  +  Sv  (4) 

where  /z  =  S\ii\  +  (1  -  S\)/Jbg  (5) 

For  a  mixture  of  liquid  and  gas  phases,  the  local  mass  average 
velocity  v  is  defined  as  [33,34]: 

v_  SiAVi  +  SgPgVg  (6) 

Si  Pi  +  SgPg 

The  source  term  for  the  momentum  equation  used  in  the  model 
describes  the  flow  of  the  fluid  through  a  porous  media  by  using 
Darcy’s  drag  force.  The  gravity  and  surface  tension  forces  were  also 
considered  in  the  momentum  equation  source  term.  The  source 
terms  for  different  regions  of  the  fuel  cell  are  given  in  Table  2. 

2.2.3.  Volume  fraction  equation 

The  tracking  of  the  interface  between  the  phases  was  accom¬ 
plished  by  the  solution  of  a  continuity  equation  for  the  volume 
fraction  of  one  (or  more)  of  the  phases.  The  motion  of  the  interface 
between  two  immiscible  liquids  (namely,  gas  mixture  and  liquid 
water)  of  different  density  and  viscosity  in  the  VOF  method  was 
defined  by  volume  fraction  of  liquid  water  s\  and  volume  fraction 
of  gaseous  phase  sg  [33]. 

For  volume  fraction  of  liquid  water: 


^(esiPi)  +  v(siPiv,)  =  Ss 


(7) 


2.2.4.  Energy  conservation  equation 

In  this  model,  the  energy  balance  in  terms  of  temperature 
change  was  also  considered.  In  the  multi-phase  model,  the  energy 
equation  is  also  shared  among  the  phases  [33,35]: 


(/°cp)eff§-  +  (pcp)eff(vvr)  =  V(keffVT)  +  Sr 
where  T  is  the  temperature  (I<)  that  is  defined  as 


T  = 


siPih  +  (1  -  si)Pg7g 

SiPi+(l  —  si)pg 


(8) 


(9) 


To  specify  these  parameters  in  porous  media,  the  effective  prop¬ 
erties  were  determined: 


(Pcp)e  ff  —  £Pfcp,f  +  (1  —  £)PsCp,s 

1<eff  =  £k$  +  (1  —  &)kf 


(10) 

(11) 


The  source  term  of  energy  equation  for  two-phase  model  may 
consist  of  heat  from  electrochemical  reactions,  heat  of  formation 


Table  2 

The  source  terms  of  governing  equations 


Governing  equation 

Volumetric  source  terms  and  location  of 
application 

Conservation  of  mass 

For  all  parts:  Sm  =  0 

Volume  fraction 

For  all  parts:  Ss  =  rw,  rw  = 
cr  max  [(1  -S|)PwYTPsatMH2o,  -Sja] 

Conservation  of  momentum 

For  gas  channels:  Sv  =  0 

For  GDLs  and  void  of  catalyst  layers: 

Sv  =  pg  ^v  +  XK,^') 

For  membrane: 

Sv  =  pg  te2v  +  XK  +  *  cfnfFV</>m 

Conservation  of  energy 

For  current  collectors:  St  = 

For  gas  flow  channel:  Sj  =  0 

For  GDL:  ST  =  +  rw/iL 

For  membrane:  Sj  =  ^ 

For  catalyst  layer: 

/  \ 

St  =  pRan,cat  +  f 2  f  J  +  h/v^L 

Hydrogen  transport 

Mh 

For  anode  catalyst  layer:  Sh2  =  -  -jfrR an 

Oxygen  transport 

For  cathode  catalyst  layer:  So2  =  -^jrRcat 

Water  vapor  transport 

For  anode  catalyst  layer: 

/  ndMH  0\ 

Sh20  —  [  p  1  ^an  G/v 

For  cathode  catalyst  layer: 

Sh20  =M^°  Rea, +  ("d>°)  Scat -rw 

Conservation  of  Charge 

For  anode  catalyst  layer: 

Ses  —  —Ran',  $em  —  Ran 

For  cathode  catalyst  layer: 

Ses  =  ^catl  Sem  =  —Rc at 

For  other  parts:  Ses  =0;  Sem  =  0 

of  water,  Ohmic  heating  caused  by  the  Ohmic  resistance  of  solid 
phases  and  the  heat  due  to  phase  change.  For  different  zones,  the 
source  terms  are  expressed  as  in  Table  2. 


2.2.5.  Species  transport  equations 

The  model  predicts  the  local  mass  fraction  of  each  species,  Y2-, 
through  the  solution  of  a  convection-diffusion  equation  for  the 
ith  species.  The  species  transport  equations  are  generally  in  the 
following  form: 

jfzpY,)  +  V  ■  (spvYj)  =  Du  mV2(pY,)  +  Si  (12) 


where  Di  m  is  the  diffusion  coefficient  for  species  i  in  the  mixture 
[34], 


Di,m  =  elb(l-Si)rsD° 


(13) 


where  D°m  is  the  diffusion  coefficient  for  species  i  in  the  mix¬ 
ture  at  reference  temperature  and  pressure.  The  source  term  in  the 
transport  equation  is  shown  in  Table  2. 

During  the  operation,  the  H+  move  from  the  anode  to  the  cath¬ 
ode  and  also  pull  water  molecules  with  them,  this  is  known  as 
the  electro-osmotic  drag  effect.  Physically,  the  water  transport  rate 
through  the  membrane  from  anode  to  cathode  by  electro-osmotic 
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drag  is  computed  as 

^d^H20  , 


mH20  = 


-Feat 


(14) 


where  nd  is  the  drag  coefficient  and  is  proposed  by  Springer  et  al. 
[36]  for  Nafion  membrane: 


(15) 


where  X  is  the  water  content  inside  the  polymer  membrane  [36]. 


2.2.6.  Conservation  of  charge 

The  current  transport  of  electrons  through  the  solid  phase  and 
ions  through  the  membrane  phase  was  represented  by  the  follow¬ 
ing  equations  [33,35]: 


2.3.  Water  transport  and  its  effect  on  the  properties  of  porous 
media 

Water  is  formed  in  the  cathode  catalyst  layer  by  electrochem¬ 
ical  reactions— an  amount  of  oxygen  is  consumed  and  an  amount 
of  water  is  produced.  Due  to  the  proton  movement  from  the  anode 
to  the  cathode  through  the  membrane,  water  molecules  are  pulled 
with  the  protons  by  a  force  called  electro-osmotic  drag  [30].  Addi¬ 
tionally,  water  may  diffuse  through  the  membrane  due  to  the 
concentration  differences. 

In  the  present  study,  the  membrane  water  diffusivity  and  mem¬ 
brane  ionic  conductivity  were  calculated  by  the  expressions  given 
by  Springer  et  al.  [36].  It  is  almost  a  standard  model  that  has  been 
well  explained  in  the  recent  textbook  [37]. 


V  •  (crsV0s)  —  Ses 


(16) 


2.4.  Volume-of-fluid  (VOF)  model 


V  •  (<7m  V0m)  —  Sem 


(17) 


The  volumetric  source  terms  Ses  and  Sem  are  defined  as  volu¬ 
metric  transfer  currents.  The  source  term  in  the  equation  of  charge 
is  shown  in  Table  2. 

The  source  terms  representing  transfer  current  were  calculated 
by  using  Butler-Volmer  equation  [33]: 


The  relation  between  species  concentration  and  species  mass 
fraction  is 


(20) 


The  volumetric  transfer  current  R  is  driven  by  the  activation 
overpotential  p,  which  is  the  potential  difference  between  solid  and 
membrane  phases  [33,35]: 


h  —  0s  -  0m  -  Kef 


(21) 


The  reference  potential  of  the  electrode  Kef  is  0  on  the  anode 
side  and  is  equal  to  the  open  circuit  voltage  on  the  cathode  side: 


Kef  =  0,  on  anode  side;  Kef  =  K>c,  on  cathode  side  (22) 

The  cell  potential  is  then  the  difference  between  the  cathode 
and  anode  solid  phase  at  the  terminal  collectors  (two  ends  of  the 
cell  collectors  which  are  connected  to  the  external  electric  circuit 
from  the  both  electrodes). 


Kell  —  0s,  cat  —  0s,  an 


(23) 


The  VOF  technique  was  implemented  in  the  channels  and  porous 
media  (including  GDLs  and  catalyst  layers)  [33].  Interface  between 
gas  and  liquid  (two-phase  flow)  is  tracked  by  the  volume  fraction  of 
liquid  water  in  the  computational  cell  volume.  In  the  VOF  approach, 
the  source  term  of  continuity  and  momentum  equations  used  in 
porous  media  include  the  effects  of  surface  tension  and  wall  adhe¬ 
sion  and  capillary  water  transport  phenomenon,  therefore,  was  also 
considered. 


2.4.1  Geometric  reconstruction  scheme 

The  geometric  reconstruction  PLIC  scheme  (piecewise  linear 
interface  construction)  was  employed  because  of  its  accuracy  and 
applicability  for  general  unstructured  meshes,  compared  to  other 
methods  such  as  the  donor-acceptor,  Euler  explicit,  and  implicit 
schemes.  A  VOF  geometric  reconstruction  scheme  is  divided  into 
two  parts:  a  reconstruction  step  and  a  propagation  step.  Details 
can  be  found  in  Refs.  [33,38]. 


2.4.2.  Implementation  of  surface  tension 

The  addition  of  surface  tension  to  the  VOF  method  is  modeled 
by  a  source  term  in  the  momentum  equation.  The  pressure  drop 
across  the  surface  depends  upon  the  surface  tension  coefficient  x 
[33,38]: 

a|,-*(s;  +  i5)  (26) 


where  Ri  and  R2  are  the  two  radii,  in  orthogonal  directions,  to 
measure  the  surface  curvature. 

The  surface  tension  can  be  written  in  terms  of  the  pressure  jump 
across  the  interface,  which  is  expressed  as  a  volumetric  force  F 
added  to  the  momentum  equation: 


F  =  xk 


pv*i 

(p\  +  pg)/2 


(27) 


In  order  to  satisfy  the  conservation  of  charge,  the  total  current 
of  either  electrons  or  protons  coming  out  from  the  anode  catalyst 
layer  must  be  equal  to  the  total  current  coming  into  the  cathode 
catalyst  layer  and  must  be  equal  to  the  total  current  caused  by  the 
proton  movement  through  the  membrane  [35]: 


where  x  denotes  the  surface  tension  coefficient  and  k  is  the  surface 
curvature.  The  source  terms  for  different  regions  of  the  fuel  cell  are 
given  in  Table  2. 

The  surface  curvature  k  in  Eq.  (27)  can  be  defined  in  terms  of  the 
divergence  of  the  normal  unit  vector  of  the  interface  h: 


van  vcat  (24) 


The  average  current  density  is  the  total  current  generated  in  the 
fuel  cell  divided  by  the  geometric  area: 


(25) 


k  =  V  •  h  =  V  •  (hwcos^w  +  twSin^w)  (28) 

where  h  is  the  unit  vector  normal  to  the  interface  between  two 
phases  near  the  walls,  hw  is  the  unit  vector  normal  to  the  walls,  fw 
is  the  unit  vector  tangential  to  the  walls,  and  0W  is  the  static  contact 
angle  at  the  walls  shown  in  Fig.  2.  For  the  electrode  surfaces  with 
different  wettabilities,  different  static  contact  angles  were  assigned, 
and  different  contact  angles  could  result  in  different  surface  ten- 
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air 


interface 


Qq2,  in 


=  (Oq2  Mq2 


^ave^surf 

4Fpo2Ain 


(32) 


2. 6.2.  Ou tlet  of  flow  channels 

The  fully  developed  condition  was  assumed  to  be  applied  for 
velocity  field  and  species  concentrations. 


Fig.  2.  Static  contact  angle  at  the  walls. 

sions  (F),  thus  influencing  on  water  transport.  The  values  of  surface 
tension  and  contact  angle  used  in  this  model  are  listed  in  Table  1. 

2.5.  Mass  transfer 


2.6.3.  The  terminal  collectors 

There  were  two  terminal  collectors  in  which  the  anode  and  cath¬ 
ode  collectors  were  connected  to  the  external  electric  circuit.  It  was 
necessary  to  determine  the  boundary  conditions  for  phase  poten¬ 
tials  on  these  interfaces. 


Liquid  water  is  considered  to  appear  in  the  fuel  cell  when  the 
water  vapor  pressure  reaches  its  saturated  value  at  the  cell  oper¬ 
ating  temperature.  Contrarily,  liquid  water  is  evaporated  when  the 
water  vapor  pressure  is  smaller  than  its  saturated  value.  Water  con¬ 
densation  and  evaporation  phase  change  is  also  an  important  factor 
to  determine  the  presence  of  liquid  water  in  multi-phase  model.  In 
the  present  study,  the  variation  of  water  mass  flow  rate  was  due  to 
condensation  (or  evaporation)  as  described  here  [33]: 


fn  =  rw 


Mh2o(A/v  —  Psat) 

Kf 


(29) 


2.6.4.  External  boundaries 

External  boundaries  are  defined  as  all  outside  surfaces  of  PEMFC 
except  the  terminals  (the  outside  surfaces  of  current  collectors 
where  the  current  entry  to  or  departure  from).  Zero-current-flux 
condition  was  applied  for  the  external  boundaries  due  to  no  cur¬ 
rents  coming  or  leaving. 

For  the  sake  of  simplicity,  the  cooling  channel  was  not  included 
in  the  present  model.  However,  it  is  necessary  to  control  the  cell 
temperature.  Consequently,  the  heat  transfer  was  assumed  to  take 
place  between  the  all  outside  surfaces  and  the  ambient  environ¬ 
ment  by  forced  convection  as  shown: 


2.6.  Boundary  conditions 


Qconvection  —  hAsi^s  T0C)(W) 


(33) 


To  close  the  equation  systems  including  conservation  equa¬ 
tions  of  continuity,  momentum,  energy,  species,  and  charge  with 
unknowns:  u,  v,  w,  p,  T,  Y,  and  0,  the  boundary  conditions  are 
required.  The  summary  of  boundary  conditions  is  listed  in  Table  3. 
Below  a  brief  description  is  provided. 


where  h  is  the  convection  heat  transfer  coefficient  (W  rrr2  K-1 ),  As 
is  heat  transfer  surface  area  (m2),  Ts  and  are  the  temperatures 
of  the  surfaces  and  the  free  stream,  respectively. 

Noted,  heat  is  produced  when  the  fuel  cell  operates.  For  the  case 
in  which  water  finally  ends  in  vapor  form,  this  heat  generation  rate 
is  defined  as  [39]: 


2.6.1.  Inlet  of  flow  channels 

Inlet  velocities,  fuel  and  oxidant  temperatures  and  mass  concen¬ 
tration  of  species  were  set  as  given  parameters  listed  in  Table  3.  Inlet 
velocities  can  be  preliminarily  calculated  by  the  inlet  flow  rates 
based  on  the  average  current  density  (/avg)  and  excess  coefficient  co 
that  is  defined  as 


0)h2 


^H2,  supply 
^H2 ,  consumption 


^02 


m02,  supply 
™02 ,  consumption 


(30) 


Generated  =/(1.25-Vcell)(W)  (34) 

In  order  to  ensure  that  the  forced  convection  dissipates  the  heat 
converted  from  electricity,  the  following  condition  must  be  satis¬ 
fied: 


Qconvection  —  Qgenerated  (35) 

then  the  convective  coefficient  h  was  chosen  as 


Then,  inlet  velocity  expression  is  given  by  the  following  equations: 


u  A/r  ^aveAsurf 

UH2,,n  H2MH2  2FpWlAm 

(31) 

Table  3 

Boundary  conditions 

Locations  of  application 

Boundary  conditions 

Inlet  of  the  anode  flow  channel 

u  =  ^H2  =  ^H2,in!  ^H20 

=  ^H20,an,in»  Fan, in  =  F-l2,in 

Inlet  of  the  cathode  flow  channel 

u  =  ^cat,in!  Yd2  —  Yo2,in»  ^H20 

=  ^H20,cat,in*  F:at,in  =  Fair/02,in 

Outlet  of  the  anode  flow  channel 

3u0ut,an  n.  9Yh2  p.  aYH20,an 

dx  ~  dx  ~  U’  8x 

r\.  3ran,out  r\ 

~  u’  dx  - 

Outlet  of  the  cathode  flow  channel 

9liout,cat  p.  9Y02  p.  3yH20,cat 

dx  ’  dx  — u’  dx 

r\.  9Tcat,out  n 
“  U’  dx  ~U 

The  anode  terminal 

II 

o 

II 

O 

The  cathode  terminal 

II 

2. 

&\f 

II 

O 

External  boundaries 

90s  p .  90s  p .  90m  p .  90m 
dx  -u’  dz  -u'  dx  -u’  dz 

=  0 

7(1.25 -ycell) 
A^s-Too) 


(36) 


2.7.  Solution  procedure 

The  above-coupled  set  of  governing  equations  and  relative  equa¬ 
tions  were  implemented  into  FLUENT®  6.2  by  developing  our  own 
user-defined-functions  (UDFs)  based  on  the  general  FLUENT®  pack¬ 
age  that  does  not  include  the  PEMFC  module  developed  by  Fluent 
Inc.  The  developed  own  UDFs  are  written  in  C  language  with 
about  2000  statements.  There  were  totally  279,500  grid  cells  in 
the  computation  domain  used  to  simulate  the  physical  and  elec¬ 
trochemical  phenomena  in  the  fuel  cell.  The  solution  procedure  for 
pressure-velocity  coupling  was  based  on  PISO  algorithm  [40]. 

In  order  to  save  the  calculation  time,  an  unsteady  single-phase 
model  was  simulated  from  the  time  t=0s  to  0.5  s.  At  the  time 
t= 0.5  s,  an  initialization  of  a  series  of  liquid  water  droplets  was 
initially  suspended  on  the  cathode  channel  at  the  time  t= 0.50  s 
and  a  general  model  of  PEMFC  was  called  to  further  investigate 
two-phase  flow  behavior,  especially  liquid  water  behavior  across 
porous  media,  together  with  electrochemical  reaction,  heat  and 
mass  transfer. 
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Fig.  3.  A  series  of  water  droplets  was  initially  suspended  into  the  cathode  channel. 


3.  Results  and  discussion 

Fig.  3  shows  the  arrangement  of  the  initial  liquid  water  distri¬ 
bution.  There  were  15  equidistant  spherical  droplets  (r=0.4mm) 
freely  suspended  along  the  cathode  channel.  The  cell  voltage  was 
constantly  set  to  be  0.5  V.  The  operating  condition  and  fuel  cell 
parameters  are  listed  in  Table  1. 

3.1.  Water  distribution  inside  the  fuel  cell— comparison  of 
numerical  simulation  and  experimental  visualization 

Before  we  have  comparisons,  let  us  have  a  quick  look  at  the 
experimental  investigations  in  the  available  literature.  Experimen¬ 
tal  studies  to  probe  detailed  liquid  water  transport  from  the  GDL 
into  the  gas  flow  channels  have  been  performed  by  Yang  et  al.  [41  ] 
and  Zhang  et  al.  [42].  Three  transparent  PEMFCs  with  different  flow 
fields,  including  parallel,  interdigitated  and  cascade,  were  intro¬ 
duced  by  Liu  et  al.  [43].  As  presented  in  Liu’s  paper,  the  effects 
of  flow  field,  cell  temperature,  cathode  gas  flow  rate  and  opera¬ 
tion  time  on  water  build-up  and  cell  performance  were  presented, 
respectively.  Elowever,  only  the  water  flooding  and  two-phase  flow 
in  the  transparent  channel  is  visualized.  Recently,  there  are  some 
innovative  visualization  methods  that  were  applied  to  investigate 


Fig.  4.  (a)  Water  concentration  distributions  in  the  porous  media  on  the  cross- 
sectional  planes  (Y-Z  plane)  at  X= 0.01  m;  (b)  water  profile  in  two  2.25  cm2  active 
area  cells  with  50%  RH,  1.1  stoich  and  172  kPa  at  the  anode  and  2.2  stoich  and  172  kPa 
at  the  cathode  at  1.4  A  cm-2  and  100%  RH  at  the  cathode  inlet  (Davey  et  al.  [46]). 


in  situ  the  production  and  distribution  of  water  in  a  PEMFC.  Fein- 
del  et  al.  [44]  developed  a  nuclear  magnetic  resonance  (NMR) 
microscopy  system  to  investigate  the  distribution  of  water  through¬ 
out  an  operating  H2/02  PEMFC.  The  areas  investigated  include 
the  flow  fields,  current  collectors,  membrane  electrode  assem¬ 
bly  (MEA),  and  the  membrane  surrounding  the  catalysts.  Another 
method  is  neutron  radiography  and  tomography  that  has  been  suc¬ 
cessfully  applied  for  investigation  of  liquid  water  evolution  and 
transport  in  low  temperature  PEMFCs  by  Manke  et  al.  [45]  and 
Davey  et  al.  [46]. 

Now  let  us  make  the  comparisons.  The  water  distribution  in 
the  channel  and  porous  media  of  the  fuel  cell  could  be  obtained 
from  the  numerical  model.  To  validate  the  numerical  results  pre¬ 
sented  in  this  study,  qualitative  comparisons  were  made  by  using 
some  available  experimental  visualization  data.  Fig.  4a  presents  the 
computed  water  distribution  in  a  Y-Z  cross-sectional  plane  of  this 
fuel  cell  model.  Fig.  4b  depicts  the  water  distribution  captured  by 
using  Neutron  Radiography  method  [46].  Note  that  in  Fig.  4b,  there 
was  more  water  (black  =  no  water;  blue/white  =  more  water)  accu¬ 
mulated  near  the  outlets  than  near  the  inlets.  The  water  near  the 


Fig.  5.  (a)  Volume  fraction  distribution  of  liquid  water  in  the  middle-plane  of 
the  cathode  channel  at  t= 0.502  s;  (b)  neutron  tomogram  showing  the  three- 
dimensional  water  distribution  in  the  cathode  channel  of  a  serpentine  PEMFC 
(Manke  et  al.  [45  ] ):  (b-a)  right  before  shut  up,  (b-b)  several  hours  after  a  tomography 
was  performed,  (b-c)  quotient  of  image  (b-a)  and  (b-b),  (b-d)  neutron  tomogram 
showing  the  three-dimensional  water  distribution. 
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Fig.  6.  The  motion  and  deformation  of  droplets  vs.  time  on  the  cathode  channel:  (a)  t= 0.500  s;  (b)  t= 0.50025  s;  (c)  t  =  0.5005  s;  (d)  t  =  0.50075  s;  (e)  the  velocity  distribution 
and  the  volume-of-fluid  of  liquid  water  in  theX-Z  plane  at  Y=  0.002  m  at  t= 0.50075  s;  (f)  the  pressure  distribution  and  the  volume-of-fluid  of  liquid  water  in  theX-Z  plane  at 
Y  =  0.002  m  at  t  =  0.50075  s;  (g)  the  water  and  velocity  distribution  in  the  different  cross-sections  along  the  main  flow  direction  at  t= 0.50075  s,  corresponding  to  (a)X=  0.001  m, 
(b)X=  0.006  m,  (c)X=  0.014  m  and  (d)X= 0.019  m,  respectively. 
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(9)  Reference  vector  5m/s 
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Fig.  6.  ( Continued ). 


outlet  was  more  concentrated  in  the  cathode  GDL,  especially  the 
GDL  above  the  land  area.  This  is  consistent  with  the  generation  of 
water  at  the  cathode  and  its  removal  by  the  gas  flow  in  the  channels. 
The  water  generated  above  the  land  area  needs  to  diffuse  out  later¬ 
ally  to  the  channel  before  it  can  be  removed  through  the  flow  fields. 
The  numerical  result  in  Fig.  4a  qualitatively  shows  good  agreement 
with  the  experimental  data  in  Fig.  4b  made  by  Davey  et  al.  [46].  In 
addition,  more  liquid  water  concentrated  in  the  turning  area  of  ser¬ 
pentine  channel  fuel  cell  is  also  presented  by  the  computed  model 
in  this  study  shown  in  Fig.  5a  and  observed  in  Fig.  5b  that  was  cap¬ 
tured  by  using  neutron  tomography  method  in  the  experiments  of 
Manke  et  al.  [45]. 

In  summary,  by  comparing  the  numerical  results  with  avail¬ 
able  experimental  results,  we  can  conclude  that  the  present  general 
model  is  qualitatively  in  good  agreement  with  the  available  experi¬ 
mental  results.  For  this  research,  qualitative  comparisons  could  not 
be  conducted  at  the  present  time  due  to  the  sizes  of  the  fuel  cell  in 
the  cited  experimental  research  are  not  known.  The  present  authors 
are  currently  developing  a  bench  experimental  set-up  to  further 
validate  the  general  model  presented  here  and  detailed  results  will 
be  published  in  future. 

3.2.  The  motion  and  deformation  of  liquid  water  inside  the 
cathode  channel 

Fig.  6  shows  the  motion  and  deformation  of  droplets  versus  time 
on  the  cathode  channel  of  the  PEMFC.  At  t  =  0.5  s,  the  droplet  was 
placed  at  the  channel  inlet  and  was  in  its  initial  spherical  shape 
(Fig.  6a).  As  time  progressed,  droplet  deformation  (elongation)  in 
the  flow  direction,  which  is  due  to  a  combination  of  droplet  sur¬ 
face  tension  and  shear  stress  from  surrounding  airflow,  could  be 
observed.  In  other  words,  the  enlargement  of  droplet  surface  area 
increased  the  surface  tension  such  that  the  force  balance  on  the 
droplet  could  be  achieved.  Fig.  6b  indicates  that,  after  hitting  the 
turning  surface,  these  deformed  droplets  had  the  tendency  of  frag¬ 
menting  and  then  entering  the  central  airflow  due  to  the  dragging 
effect  from  the  turning  flow.  Obviously,  in  the  turning  area,  the 
shear  force  from  the  airflow  is  much  stronger  than  the  droplet 
surface  tension,  thus  together  with  the  inertia  force,  resulting  in 
significant  droplet  deformation  and  complex  liquid  water  distri¬ 
bution.  As  the  time  progressed  after  milliseconds,  the  deformed 


droplets  that  enter  the  turning  area  from  the  straight  channels  hit 
the  turning-wall  surface.  As  a  result,  the  deformed  droplets  were 
further  elongated  and  the  strong  “turning”  gas  flow  breaks  these 
droplets  into  even  smaller  ones  sticking  on  the  turning-wall  sur¬ 
face.  Due  to  the  effect  of  wall  adhesion  and  surface  tension,  the 
smaller  water  droplets  sticking  on  the  turning-wall  would  slowly 
move  forward  while  the  other  water  droplets  from  the  channel  is 
continuously  coming  to  the  turning,  and  further  broken  into  the 
other  smaller  ones.  Such  these  water  droplets  would  coalesce,  then 
expanding  and  deforming  into  “water  bands”  due  to  the  shear  stress 
from  “turning”  gas  flow  and  pressure  drop,  as  shown  in  Fig.  6c  and 
d.  It  is  noted  that  since  the  coalescence  of  water  droplets  is  formed 
in  the  turning  area  resulting  in  a  high  concentration  of  liquid  water 
volume,  a  substantial  blockage  in  the  channel  would  take  place.  This 
blockage  induces  a  decrease  of  the  gap  between  the  water  droplets 
and  the  channel  wall.  In  addition,  the  water  droplets  hitting  the 
turning  wall  would  stick  onto  it  due  to  the  wall  adhesion.  Fig.  6e 
shows  the  velocity  distribution  and  the  volume  fraction  of  liquid 
water  in  the  X-Z  plane  at  Y=  0.002  m  at  t  =  0.50075  s.  The  velocity 
field  of  airflow  is  accelerated  through  the  gaps  that  were  formed  by 
the  blockage  of  liquid  water  in  the  turning  area,  and  then  is  decel¬ 
erated  after  leaving  the  turning  area  in  which  the  blockage  occurs. 
In  the  other  views  shown  in  Fig.  6g,  one  can  see  that  the  water 
droplets  hit  the  turning  wall  and  expanded  to  the  surroundings  of 
the  landing  points.  The  droplets  could  induce  high-pressure  zones 
in  the  channel,  as  shown  in  Fig.  6f,  especially  in  the  turning  area. 

As  the  time  progressed  to  t  =  0.501  s,  as  shown  in  Fig.  7a,  the 
water  bands  in  the  turning  walls  keep  moving  forwards  along 
the  main  flow  direction.  One  can  see  in  Fig.  7a  and  b  that,  the 
water  bands  have  a  tendency  of  elongating  without  breaking  away 
under  the  impact  of  shear  force  from  the  gas  flow.  However,  the 
water  bands  were  broken  again  into  smaller  droplets  distributed 
over  the  straight  channel  wall  since  these  bands  left  the  turn¬ 
ing  area.  Subsequently,  the  formed  small  water  drops  are  easily 
removed  and  distributed  everywhere  on  the  straight  channel  walls, 
as  it  can  be  seen  in  Fig.  7c  and  d  at  t  =  0.504  s  and  0.506  s.  Fig.  7e 
shows  that  the  airflow  is  rapidly  accelerated  at  the  turning  area 
due  to  the  water  blockage  as  earlier  mentioned.  Hence,  the  shear 
stress  from  the  airflow  in  the  straight  channel  is  significant  to 
break  the  water  bands  into  small,  discrete  liquid  droplets  rather 
than  deforming  or  elongating  such  water  bands.  This  also  results 
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Reference  vector  5m/s 


VOFofliquid water:  0.32  0.44  0.56  0.68  0.8 


Y  =  0.002  m 


Fig.  7.  The  motion  and  deformation  of  droplets  vs.  time  on  the  cathode  channel:  (a)  t=  0.501  s;  (b)  t= 0.502  s;  (c)  t= 0.504  s;  (d)  t= 0.506  s;  (e)  the  velocity  distribution  and 
the  volume-of-fluid  of  liquid  water  in  the  X-Z  plane  at  Y=  0.002  m  at  t= 0.506  s;  (f)  the  pressure  distribution  and  the  volume-of-fluid  of  liquid  water  in  the  X-Z  plane  at 
Y=  0.002  m  at  t  =  0.506  s. 
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(a) 


(b) 


0.02 


Fig.  8.  The  motion  and  deformation  of  droplets  vs.  time  on  the  cathode  channel:  (a)  t=  0.520  s;  (b)  t  =  0.522  s;  (c)  t= 0.524  s;  (d)  t=  0.526  s;  (e)  the  velocity  distribution  and 
the  volume-of- fluid  of  liquid  water  in  the  X-Z  plane  at  7=0.002  m  at  t= 0.526  s;  (f)  the  pressure  distribution  and  the  volume-of-fluid  of  liquid  water  in  the  X-Z  plane  at 
Y=  0.002  m  at  t=  0.526  s. 
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Fig.  9.  The  motion  and  deformation  of  droplets  vs.  time  on  the  cathode  channel:  (a)  t  =  0.560  s;  (b)  t= 0.562  s;  (c)  t  =  0.564  s;  (d)  t= 0.566  s;  (e)  the  velocity  distribution  and 
the  volume-of-fluid  of  liquid  water  in  the  X-Z  plane  at  Y=  0.002  m  at  t  =  0.566  s;  (f)  the  pressure  distribution  and  the  volume-of-fluid  of  liquid  water  in  the  X-Z  plane  at 
Y=  0.002  m  at  t  =  0.566  s;  (g)  the  liquid  water  in  the  Y-Z  plane  at  X= 0.001  m,  0.002  m,  0.018  m  and  0.019  m. 
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Fig.  9.  ( Continued ). 


in  a  decrease  of  pressure  drop  in  comparison  with  the  value 
at  t= 0.50075  s. 

At  t= 0.52  s  as  shown  in  Fig.  8a,  a  portion  of  water  droplets  was 
removed  from  the  cathode  channel  outlet  by  the  gas  flow.  The  small 
droplets  tend  to  scatter  everywhere  in  the  channel  rather  than  coa¬ 
lescing  or  forming  water  bands.  Even  in  the  turning  area,  it  can 
be  seen  that  the  water  bands  either  are  tiny  or  are  not  formed  at 
all.  The  wall  adhesion  of  the  small  droplets  previously  coalesced 
at  the  turning  areas  is  insufficient  to  resist  the  shear  stress  from 
the  gas  flow.  In  addition,  one  can  see  that  there  are  some  droplets 
penetrating  to  the  gas  diffusion  layer.  It  can  be  realized  that  this 
phenomenon  takes  place  mostly  in  the  turning  areas  where  the 
water  is  not  easy  to  be  removed  rather  than  in  the  straight  chan¬ 
nel  as  previously  discussed.  As  shown  in  Fig.  8b-d  that  as  the  time 
proceeds  to  t  =  0.526  s,  the  water  continues  to  be  taken  away  from 
the  channel  outlet  except  the  water  amount  remaining  in  the  gas 
diffusion  layer.  At  t  =  0.526  s,  as  shown  in  Fig.  8e,  the  velocity  field 
is  observed  to  be  “more  stable”  since  a  large  amount  of  liquid  water 
was  taken  away  and  then  the  blockage  caused  by  the  small  droplets 
less  influences  on  the  gas  flow  movement.  Consequently,  this  also 


results  in  an  evident  decrease  in  pressure  drop  distribution  which 
is  presented  in  Fig.  8f. 

Most  of  the  water  droplets  were  moved  out  of  the  channel  by  the 
gas  flow  at  t=  0.56  s.  As  the  time  proceeded  to  t  =  0.566  s,  it  could 
be  observed  that  there  were  only  few  small  droplets  remaining 
in  the  exit  channel.  Flowever,  an  amount  of  liquid  water  pene¬ 
trated  through  the  gas  diffusion  layer  seemed  to  be  retained  at 
the  same  position  in  the  layer  regardless  the  droplets  in  the  chan¬ 
nel  continued  to  be  removed  outwards.  It  can  be  obviously  seen 
that  the  water  amount  kept  in  the  diffusion  layer  is  marked  by  cir¬ 
cles  as  shown  in  Fig.  9a-d.  Noticeably,  the  marked  circles  in  Fig.  9 
evoke  that  the  penetrated  water  usually  occur  in  the  turning  area 
in  which  the  liquid  water  is  more  difficult  to  be  removed  and  have 
a  tendency  of  passing  through  the  gas  diffusion  layer  instead  of 
moving  along  the  channel  to  the  outlet.  The  liquid  water  in  the 
gas  diffusion  layer  is  also  clearly  illustrated  in  Fig.  9g  since  the 
cross-sectional  planes  along  X-direction  at  X=  0.001  m,  0.002  m, 
0.018  m  and  0.019  m  are  extracted.  As  time  progresses,  the  more 
liquid  water  gradually  move  out  of  the  channel,  the  less  the  liquid 
water  blocks  the  flow  field  in  the  channel.  Therefore,  the  block- 
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Fig.  10.  (a)  The  distribution  of  volume  fraction  of  liquid  water  in  the  cross-sectional  planes  (Y-Z  plane)  at  X=  0.0015  m,  0.010  m  and  0.0185  m  (Case  2);  (b)  distribution  of 
volume  fraction  of  liquid  water  in  the  middle-plane  of  the  cathode  flow  channel  (Y=  0.0015  m)  and  catalyst  layer  (Y=  0.002305  m)  (Case  2). 
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Fig.  11.  Velocity  vectors  in  the  channel:  on  the  middle-plane  of  cathode  channel  (Y=  0.0015  m)  and  the  middle-plane  of  anode  channel  (Y=  0.00317  m):  (a)  Case  1  and  (b) 
Case  2;  on  the  cross-sectional  planes  (Y-Z  plane)  atX= 0.0015  m,  0.010  m  and  0.0185  m:  (c)  Case  1  and  (d)  Case  2. 
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Fig.  11.  ( Continued ). 
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Fig.  12.  Velocity  vectors  in  the  porous  media  on  the  cross-sectional  planes  ( Y-Z  plane)  at  X=  0.0015  m,  0.010  m  and  0.0185  m:  (a)  Case  1  and  (b)  Case  2. 
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age  insufficiently  influences  on  the  velocity  and  pressure  fields  in 
the  fuel  cell  channel.  It  can  be  seen  in  Fig.  9e  and  f,  the  velocity 
and  pressure  fields  at  t= 0.566  s  are  more  stable,  have  no  sudden 
changes  and  have  a  smooth  distribution  along  the  channel,  similar 
to  those  when  compared  with  the  case  without  liquid  water  in  the 
channel. 


3.3.  Effects  of  liquid  water  on  the  fuel  cell  characteristics 

In  order  to  investigate  the  influences  of  liquid  water  on  the  cell 
characteristics  and  operating  condition,  two  cases  (Cases  1  and  2) 
are  compared  in  this  study.  The  first  case  would  be  considered 
at  the  time  instant  t= 0.53  s  when  the  liquid  droplets  were  not 
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Fig.  13.  Velocity  vectors  and  pressure  distributions  on  the  middle-plane  of  catalyst  layers  (0.002305  m):  (a)  Case  1  and  (b)  Case  2. 
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introduced  in  the  cathode  channel  and  therefore  the  fuel  cell  oper¬ 
ates  in  the  condition  without  liquid  water.  The  second  case  would 
be  presented  at  the  time  instant  t  =  0.53  s  since  a  series  of  water 
droplets  had  been  introduced  at  t  =  0.5  s  into  the  cathode  channel 
and  evolved  for  0.03  s  that  can  be  used  to  show  the  liquid  water 


effects.  In  Case  2,  fuel  cell  characteristics  such  as  velocity,  pressure, 
temperature,  species  concentrations  and  the  fuel  cell  performance, 
etc.  are  considered  as  taking  the  appearance  of  water  droplets  into 
account.  Hence,  the  numerical  model  at  the  first  case  was  consid¬ 
ered  as  single-phase  model  and  the  second  case  was  considered 
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Fig.  14.  Velocity  vectors  and  pressure  distributions  on  the  middle-plane  of  anode  channel  (7=0.00317  m)  and  catalyst  layers  (0.002365  m):  (a)  Case  1  and  (b)  Case  2. 
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as  the  general  model  that  includes  two-phase  interface  tracking 
together  with  all  other  detailed  sub-models. 

3.3 A.  Volume  fraction  of  liquid  water 

In  Case  2,  Fig.  10a  shows  the  distribution  of  volume  frac¬ 
tion  of  liquid  water  in  the  cross-sectional  planes  (Y-Z  plane)  at 


X=  0.0015  m,  0.010  m  and  0.0185  m  and  Fig.  10b  depicts  in  the 
middle-plane  of  cathode  flow  channel  and  catalyst  layer,  respec¬ 
tively.  Note  that  the  liquid  water  was  not  taken  into  account  in 
single-phase  flow  that  was  applied  for  the  fuel  cell  prior  to  t  =  0.50  s. 
Therefore,  the  presence  of  volume  fraction  of  liquid  (liquid  water 
distribution)  is  not  applicable  to  Case  1.  It  could  be  observed  that 
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Fig.  15.  Mass  fraction  distributions  of  oxygen  on  the  middle-plane  of  cathode  channel  (Y=  0.0015  m)  and  catalyst  layers  (0.002305  m):  (a)  Case  1  and  (b)  Case  2. 
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the  liquid  droplets  distributed  in  the  fuel  cell  at  a  specified  time 
were  presented  via  volume  fraction  of  water.  To  investigate  the 
motion  and  interaction  of  liquid  droplets  inside  the  fuel  cell  and 
the  influence  of  liquid  phase  on  the  fuel  cell  operating  condition, 
the  results  of  this  simulation  would  be  analyzed  as  the  time  pro¬ 


gressed  after  0.03  s  at  t= 0.53  s.  The  change  of  the  water  occupation 
fraction  in  the  channel  and  liquid  water  transport  from  the  cath¬ 
ode  channel  to  the  porous  media  were  studied.  By  tracking  the 
presence  of  volume  fraction  of  water,  one  can  determine  where 
the  liquid  water  occurs  and  the  effects  of  presence  of  liquid  water 
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Fig.  16.  Mass  fraction  distributions  of  hydrogen  in  the  middle-planes  of  anode  channel  (7=0.00317  m)  and  catalyst  layers  (0.002365  m):  (a)  Case  1  and  (b)  Case  2. 
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on  the  physical  parameters  and  characteristics  of  the  fuel  cell. 
Fig.  10a  and  b  shows  that  most  of  the  liquid  water  concentrated 
in  the  channel  and  moved  through  porous  media  is  in  the  turning 
area  of  the  serpentine-shaped  flow  channel.  As  a  result,  this  would 
especially  influence  the  other  flow  parameters  of  the  turning  area 
also. 


3.3.2.  Liquid  water  effects  on  velocity  field  and  pressure 
distribution 

The  velocities  at  the  cathode  and  anode  inlets  are  7.0  m  s-1  and 
4.5 ms-1,  respectively.  Fig.  11a  and  b  shows  that  the  gases  (air 
and  water  vapor  at  cathode;  and  hydrogen  and  water  vapor  at 
anode)  are  supplied  from  the  inlet,  moving  along  each  branch  of  the 
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Fig.  17.  Mass  fraction  distributions  of  water  vapor  on  the  middle-plane  of  cathode  channel  (Y=  0.0015  m)  and  catalyst  layers  (0.002305  m):  (a)  Case  1  and  (b)  Case  2. 
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serpentine-shaped  channel  for  Case  1  and  Case  2,  respectively,  dra¬ 
matically  varied  their  direction  when  approaching  the  turn,  then 
go  to  the  next  branches.  In  Case  1  (Fig.  11c),  the  velocity  vectors  on 
Y-Z  planes  on  cross-sectional  planes  of  channel  branches  shown  in 
Fig.  11c  seem  to  have  similar  variation  to  each  other.  In  other  words, 
these  velocity  distributions  on  the  same  cross-section  through  all 
channel  branches  are  the  similar.  This  is  because  that  liquid  water 
was  not  considered  in  Case  1  and  therefore  there  are  no  liquid 
droplets  effects  on  the  flow  field.  In  Case  2  (Fig.  lid),  the  velocity 
and  pressure  distributions  are  quite  different.  Due  to  the  presence 
and  blockage  of  water  droplets  over  the  channel  as  previously  men¬ 
tioned,  the  flow  field  was  not  similar  to  each  other  depending  on 
where  the  droplets  located.  It  is  noticeable  that  the  pressure  drop 
in  Case  2  was  observed  to  be  higher  than  that  in  Case  1,  especially 
around  the  location  of  the  droplets,  the  pressure  drop  reaches  high¬ 
est  values  due  to  the  momentum  resistance  caused  by  liquid  water 
in  airflow. 

The  velocity  in  the  gas  channels  on  both  sides  shown  in  Fig.  10a 
and  b  were  very  high  (the  order  of  magnitude  is  1 )  by  comparing 
to  the  low  velocity  in  the  porous  media  shown  in  Fig.  12a  and  b 
(the  order  of  magnitude  is  10-2)  due  to  the  flow  resistance  in  the 
porous  media  (including  GDLs  and  catalyst  layers).  Although  the 
porous  media  is  solid,  there  still  exists  the  flow  field  due  to  pres¬ 
sure  gradient  distribution  in  the  void  of  such  regions  where  reactant 
gases  and  water  can  move  through.  It  can  be  seen  in  Fig.  12a  and  b 
for  both  cases  that  the  flow  in  Y-Z  planes  had  a  tendency  of  moving 
through  the  porous  media  from  one  channel  branch  to  the  next, 
especially  in  the  land  area  (the  porous  zones  underneath  the  col¬ 
lector  ribs).  The  difference  of  pressure  on  the  same  cross-section 
in  all  branches  is  the  driving  force  to  such  flow  motion.  Similarly, 
it  can  be  clearly  seen  in  Fig.  12b  how  the  blockage  of  liquid  water 
in  the  porous  media  influenced  the  velocity  fields.  On  the  contrary, 
for  the  case  where  liquid  water  was  absent  in  the  fuel  cell  as  shown 
in  Fig.  12a,  it  is  most  likely  that  the  flow  in  Y-Z  planes  in  the  porous 
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media  had  no  change  on  the  same  cross-section  for  every  channel 
branch. 

Fig.  13a  and  b  shows  velocity  vectors  and  pressure  distribu¬ 
tions  on  the  middle-plane  of  catalyst  layers  (0.002305  m)  for  Cases 
1  and  2,  respectively.  Now,  the  flow  field  on  X-Z  planes  could  be 
observed.  In  both  cases,  there  remained  alternate  velocity  distribu¬ 
tions  having  “triangle-shaped”  on  the  porous  media  in  which  the 
magnitude  of  velocity  vectors  increased  or  decreased  along  the  X- 
direction.  This  is  the  nature  of  fluid  flow.  It  could  be  noticed  that 
the  high  velocity  vectors  were  located  in  the  region  in  the  land 
area  between  the  two  corners  of  two  sequent  branches,  also  having 
the  largest  pressure  difference.  The  flow  in  the  corner  of  a  branch, 
therefore,  tends  to  move  towards  the  corner  of  the  next  branch  at 
the  same  side  through  the  porous  media  instead  of  moving  along 
the  channel.  This  is  the  shortest  way  and  easiest  way  as  well.  The 
magnitude  of  such  velocity  vectors  decreases  along  the  channel  due 
to  the  decrease  of  pressure  difference  and  it  is  lowest  at  the  corner 
at  opposite  side  where  the  flow  freely  moves  in  the  channel  turn. 
Moreover,  the  chaos  of  velocity  distribution  caused  by  the  effect  of 
liquid  water  corresponding  to  the  location  of  droplets  in  the  porous 
media  is  more  realized  as  shown  in  Fig.  13b. 

Fig.  14a  and  b  shows  velocity  vectors  and  pressure  distributions 
on  the  middle-plane  of  anode  channel  (Y=  0.00317  m)  and  catalyst 
layers  (0.002365  m)  for  Cases  1  and  2.  Note  that  the  flow  field  in 
the  fuel  cell  may  vary  at  the  different  times.  Because  the  flow  con¬ 
ditions  at  the  anode  side  reached  the  steady  state  at  Case  1  and  the 
time  period  between  Case  1  (at  t  =  0.50  s)  and  Case  2  (at  t= 0.53  s) 
is  almost  negligible,  then  the  velocity  and  pressure  distributions 
were  almost  the  same  for  the  both  cases  due  to  the  absence  of  liquid 
water  at  the  anode  side. 

3.3.3.  Liquid  water  effects  on  concentrations 

In  the  both  cases,  it  could  be  observed  that  the  oxygen  fraction 
decreased  along  the  cathode  channel  from  the  inlet  to  the  outlet 
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Fig.  18.  The  distribution  of  local  current  density  at  the  catalyst/GDL  interfaces  at  the  cathode  and  anode:  (a)  Case  1  and  (b)  Case  2. 
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and  from  the  channel  to  the  catalyst  layer  due  to  the  oxygen  con¬ 
sumption  on  the  catalyst  layer.  Mass  transport  of  oxygen  from  the 
channel  to  the  GDL  was  mainly  driven  by  its  concentration  gradi¬ 
ent  (diffusion)  and  this  effect  was  also  partially  influenced  by  the 
velocity  component  (convection)  inside  the  porous  media.  The  flow 
in  cathode  channel  enhanced  the  gas  transport  in  the  region  under 
the  channel,  rather  than  in  the  region  under  the  rib  (the  land  area). 
As  shown  in  Fig.  15a  and  b,  high  oxygen  region  under  the  cath¬ 
ode  channel  was  observed.  On  the  other  hand,  the  land  area  had  a 
low  amount  of  oxygen.  There  are  two  reasons  accounting  for  this; 
the  first  is  that  mass  transport  of  oxygen  in  such  area  was  low  due 
to  weak  influence  of  convection  and  diffusion.  The  second  is  that 
a  large  amount  of  oxygen  was  consumed  due  to  a  high  reaction 
rate— the  chemical  reaction  had  strongly  taken  place  in  the  land 
area  which  was  directly  connected  to  the  collector  ribs  that  the 
current  density  vectors  mainly  went  through.  This  phenomenon 
will  be  discussed  in  later  section.  The  presence  of  liquid  water  sig¬ 
nificantly  influences  the  oxygen  concentration  in  the  channel  and 
porous  media  as  well.  This  is  very  important  in  water  management 
of  a  fuel  cell.  As  mentioned  in  many  studies,  liquid  water  occurring 
in  the  fuel  cell  causes  the  flooding  phenomenon  and  further  blocks 
the  oxygen  transport  through  the  GDL  and  catalyst  layer,  then  sig¬ 
nificantly  influences  the  fuel  cell  performance.  In  Case  2,  there  was 
a  noticeable  blockage  of  liquid  water  in  oxygen  concentration,  espe¬ 
cially  in  the  land  region  underneath  the  turning  area.  In  Case  1,  the 


oxygen  mass  fraction  distribution  in  the  porous  media  underneath 
the  turning  area  could  be  seen  to  have  a  uniform,  smooth  U-shaped 
as  shown  in  Fig.  15a.  Differently  in  Case  2,  the  shape  was  various, 
unsmoothed  and  arbitrary.  This  is  because  the  liquid  water  occu¬ 
pied  the  void  in  which  oxygen  was  supposed  to  occupy.  The  higher 
volume  fraction  of  liquid  water  occupies,  the  lower  the  oxygen  mass 
fraction  could  be  observed. 

Similar  to  the  oxygen  mass  fraction  distribution  in  the  cathode, 
the  hydrogen  gradually  decreased  from  the  inlet  to  the  outlet  in 
both  the  anode  flow  channel  and  porous  media  as  shown  in  Fig.  16a 
and  b.  Flowever,  the  decrease  of  hydrogen  mass  fraction  along  the 
anode  channel  was  insignificant  due  to  absence  of  liquid  water. 

For  the  both  cases,  the  water  vapor  mass  fraction  increased  along 
the  cathode  channel  by  two  sources:  by-product  water  was  gen¬ 
erated  from  the  electrochemical  reaction  at  the  cathode  and  the 
migration  of  water  from  the  anode  to  the  cathode  due  to  the  electro- 
osmotic  drag.  As  shown  in  Fig.  17  that  there  was  a  high  water  vapor 
content  remaining  in  the  land  area  at  the  cathode  as  a  result  of 
water  build-up  due  to  the  higher  reaction  rate.  Furthermore,  due  to 
the  small  velocity  field  in  the  land  area  nearby  the  edge  (shown  in 
Figs.  11  and  12),  the  accumulated  water  in  such  area  was  not  easy  to 
be  removed.  Obviously,  water  concentration  is  highest  in  the  edge 
adjacent  to  the  outlet  channel  branch  due  to  water  build-up  and 
low  convection  effect.  In  Case  2,  the  presence  of  liquid  droplets  also 
influenced  the  distribution  of  water  vapor  mass  fraction,  similar  to 
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Fig.  19.  (a)  Temperature  distribution  in  the  cross-sectional  planes  (Y-Z  plane)  atX= 0.0015  m,  0.010  m  and  0.0185  m  (Case  2);  temperature  distributions  on  the  middle-plane 
of  cathode  channel  (Y=  0.0015  m)  and  catalyst  layers  (0.002305  m):  (b)  Case  1  and  (c)  Case  2. 
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that  on  oxygen  transport.  As  can  be  seen  in  Fig.  17b,  the  higher  the 
volume  fraction  of  liquid  water  occupied,  the  larger  the  water  vapor 
mass  fraction  observed. 

3.3.4.  Liquid  water  effects  on  current  density  distributions 

Fig.  18a  and  b  shows  the  distribution  of  local  current  density  at 
the  catalyst  layer/GDL  interfaces  on  the  cathode  and  anode  sides 
for  the  both  cases,  respectively.  Note  that  the  local  current  density 
at  the  interfaces  was  generated  by  the  electron  movement  while 
the  local  current  density  at  the  mid-plane  of  the  membrane  by 
the  proton  movement  through  the  membrane.  The  current  density 
peaked  at  the  land  area  in  the  porous  media  because  the  current 
flow  tended  to  pass  through  the  land  area  of  GDLs/catalyst  layers 
which  were  in  direct  contact  with  the  collector  ribs.  Hence,  the 
current  flow  could  move  through  a  circuit  from  the  anode  to  the 
cathode  in  the  shortest  and  fastest  way.  Note  that  the  current  den¬ 
sity  is  proportional  to  the  reaction  rate.  This  also  accounts  for  the 
reason  why  reaction  rate  or  reactant  gases  consumption  was  high¬ 
est  in  the  land  area.  For  unsteady  state,  the  average  current  density 
may  vary  with  respect  to  time.  Namely,  the  average  current  den¬ 
sity  at  t  =  0.50  s  (Case  1)  was  0.624  A  cm-2  and  at  t  =  0.53  s  (Case  2) 
was  0.616  A  cm-2.  Due  to  a  small  change  of  average  current  den¬ 
sity,  it  was  not  easier  to  observe  the  differences  on  current  density 
distributions  for  both  cases. 

3.3.5.  Liquid  water  effects  on  temperature  distributions 

The  lowest  temperature  could  be  observed  at  the  flow  inlet 
where  the  flows  with  ambient  temperature  come  in,  and  the 
temperature  increases  along  the  positive  flow  direction.  The  heat 
generation  was  due  to  the  electrochemical  reaction  and  Ohmic 
heating.  Therefore,  the  highest  temperature  concentrates  in  the 
catalyst  layers  in  which  the  reactions  take  place  and  the  current 
is  generated.  In  addition,  the  temperature  distribution  through  the 
solid  materials  (including  flow  plates/collectors  and  solid  phase 
of  porous  media)  is  also  affected  by  the  heat  conduction  process. 
Depending  on  different  solid  materials,  the  thermal  conductivity 
may  change,  resulting  in  different  temperature  distributions  in  dif¬ 
ferent  solid  layers  as  shown  in  Fig.  19a.  Note  that  the  temperature 
distribution  in  the  collectors  seemed  to  be  uniform  due  to  a  high 
thermal  conductivity  of  the  collector  material.  By  controlling  the 
convective  coefficient  h  (Eq.  (36))  to  cool  the  cell,  the  average  tem¬ 
perature  was  in  the  range  of  334-340  K. 

Fig.  19b  and  c  shows  temperature  distributions  on  the  middle- 
plane  of  the  cathode  channel  (Y=  0.0015  m)  and  catalyst  layers 
(0.002305  m)  for  Cases  1  and  2,  respectively.  It  was  observed  that 
the  temperature  distribution  in  Case  1  was  higher  than  that  in  Case 
2.  The  presence  of  liquid  water  induced  a  decrease  of  temperature  at 
the  location  occupied  by  the  droplets.  Initially,  water  droplets  were 
suspended  in  the  channel  where  the  temperature  was  the  lowest  in 
the  fuel  cell  due  to  the  convection  effect.  The  temperature  of  initial 
droplets  was  set  to  be  equal  to  the  value  of  ambient  temperature 
(300 1<).  As  the  time  progressed,  these  droplets  moved  forward  as 
the  airflow  direction  in  the  channel  and  a  portion  of  liquid  water 
was  transported  through  the  porous  media.  It  is  noticeable  that 
the  temperature  in  the  porous  media  was  highest  due  to  chemical 
reactions  taking  place  on  the  catalyst  layer.  The  low  temperature  of 
liquid  water  transferred  to  and  the  heat  generated  from  the  cath¬ 
ode  catalyst  layer  resulted  in  an  uneven  temperature  distribution 
on  the  porous  media. 

3.4.  Average  current  density  with  respect  to  time 

The  average  current  density  of  the  cell  is  presented  with  respect 
to  time  as  shown  in  Fig.  20a  and  b.  For  the  time  period  from  0  s  to 
0.5  s,  the  fuel  cell  operated  without  an  addition  of  liquid  water  as 


time  (s) 

Fig.  20.  The  average  current  density  of  the  cell  vs.  time  (Case  2):  (a)  from  t=0s  to 
0.50  s  and  (b)  from  t= 0.50  s  to  0.58  s. 

shown  in  Fig.  20a.  As  the  time  progressed,  the  average  current  den¬ 
sity  slightly  increased  after  the  fuel  cell  operation  was  started.  In 
this  stage,  liquid  water  was  not  taken  into  account  and  hence  it  did 
not  influence  the  increase  of  average  current  density.  Since  a  stable 
flow  and  species  parameters  was  obtained,  it  could  be  observed  that 
the  average  current  density  reached  a  constant  value  at  6246  A  m-2 
at  t  =  0.4  s.  Fig.  20b  shows  that  when  the  water  droplets  were  added 
into  the  channel  at  t=0.5  s,  the  average  current  density  of  the  cell 
was  no  longer  constant  at  a  value  of  6246 Am-2.  In  contrast,  the 
average  current  density  tended  to  decrease.  It  could  be  realized 
that  the  presence  of  liquid  water  in  the  fuel  cell  is  a  main  reason 
to  decrease  the  current  density  of  the  cell.  As  mentioned,  water 
in  the  form  of  liquid  may  hinder  the  gas  transport,  resulting  in  a 
concentration  loss  and  therefore  cause  a  reduction  of  the  current 
density.  However,  this  decrease  was  only  maintained  until  the  time 
t= 0.508  s  at  which  the  liquid  water  was  mostly  moved  out  of  the 
channel  except  for  some  kept  in  the  porous  media.  The  absence  of 
liquid  water  led  to  a  gradual  increase  in  the  current  density  as  the 
time  progressed,  as  shown  in  Fig.  20b. 

4.  Conclusion 

In  this  study,  a  three-dimensional,  unsteady,  multi-phase,  multi- 
component  PEMFC  model  with  VOF  interface  tracking  technique 
was  solved  using  FLUENT®  by  developing  own  UDFs  to  customize 
the  governing  and  relative  equations.  Numerical  simulations  were 
conducted  for  a  single  PEFMC  with  serpentine  flow  channels.  This 
model  considered  all  the  necessary  components  consisting  of  a  unit 
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PEMFC  and  coupled  the  fluid  flow  (momentum  transport),  species 
transport,  energy  transport,  electron  and  proton  transport,  with 
volumetric  electrochemical  reactions.  The  membrane,  gas  diffusion 
layer  and  catalyst  layers  were  treated  as  volumetric  domain  with 
grid  points,  and  the  electrochemical  reactions  in  the  catalyst  layers 
as  volumetric  reactions.  Especially,  by  using  VOF  interface  tracking 
algorithm,  the  flow  behaviors  of  liquid  water  and  temperature  field 
were  also  investigated.  The  main  interpretations  and  conclusion  of 
this  study  are  summarized  as  follows: 

1.  The  motion  and  behaviors  of  liquid  water  in  a  full  3D,  unsteady 
PEMFC  was  completely  described  by  using  VOF  interface  tracking 
algorithm  combined  with  solving  the  fluid  flow,  species  trans¬ 
port,  energy  transport  and  electrochemical  governing  equations. 

2.  The  motion  of  droplets  in  the  PEMFC  model  was  illustrated  at 
different  time  steps  in  numerical  simulation.  Firstly,  the  droplets 
freely  moved  in  the  channel  then  deformed  (elongation)  in  the 
flow  direction,  which  was  due  to  a  combination  of  droplet  surface 
tension  and  shear  stress  from  surrounding  airflow.  After  hitting 
the  turning  surface,  these  deformed  droplets  had  the  tendency 
of  fragmenting  and  then  entering  the  central  airflow  due  to  the 
dragging  effect  from  the  turning  flow.  Since  the  coalescence  of 
water  droplets  was  formed  in  the  turning  area  resulting  in  a  high 
concentration  of  liquid  water  volume,  a  substantial  blockage  in 
the  channel  occurred.  After  leaving  the  turning  area,  a  significant 
shear  stress  from  the  airflow  in  the  straight  channel  broke  the 
coalescent  droplets  into  smaller,  discrete  liquid  droplets  rather 
than  deforming  or  elongating  such  water  bands.  Finally,  liquid 
droplets  were  gradually  removed  out  of  the  channel  under  a 
strong  flow  field. 

3.  The  presence  of  liquid  water  in  the  channel  significantly  influ¬ 
ences  the  flow  field  parameter.  Namely,  due  to  the  blockage  of 
water  droplets,  the  gas  flow  would  become  unevenly  distributed, 
the  high  pressure-drop  regions  would  occur  around  the  location 
of  liquid  droplets. 

4.  The  structure  of  the  channel  will  significantly  influence  the 
distribution  of  liquid  water  inside  the  cell.  As  simulated  in  a 
serpentine  channel,  the  liquid  water  is  difficult  to  be  removed 
in  the  turning  area.  Therefore,  a  portion  of  liquid  water  tends 
to  move  through  the  gas  diffusion  layer  in  the  porous  media. 
By  simulating  different  types  of  the  channels,  the  character¬ 
istics  of  liquid  water  removal  for  different  designs  could  be 
investigated. 

5.  By  quantitatively  and  qualitatively  analyzing  the  behaviors  of  liq¬ 
uid  water  and  the  average  current  densities  at  different  periods 
of  time,  it  could  be  observed  that  the  liquid  water  would  hinder 
the  gas  transport  in  the  fuel  cell,  resulting  in  a  high  concentration 
loss  and  directly  decreasing  the  current  density— in  other  words, 
severely  affect  the  fuel  cell  performance. 

6.  The  local  current  density  contours  are  also  presented,  show¬ 
ing  that  the  high  magnitude  of  the  vectors  mostly  distribute  in 
the  region  underneath  the  channel  where  high  reaction  rates 
exist. 

7.  By  solving  the  energy  equation,  the  temperature  distribution 
over  the  fuel  cell  was  described.  The  temperature  was  controlled 
at  around  333 1<  by  using  a  suitable  forced  convective  coefficient 
applied  on  the  fuel  cell.  The  model  also  showed  that  the  tem¬ 


perature  of  the  channel  and  porous  media  with  the  presence  of 
liquid  water  was  lower  than  that  without  liquid  water. 
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